Compare datasets

---

library(raster)

library(plyr)
library(dplyr)
library(reshape)

library(ggrepel)
library(ggplot2)
library(ggpmisc)
library(ggpubr)

library(viridisLite)
library(colorspace)
library(RColorBrewer)

library(sf)
library(stringr)

library(kableExtra)

memory.limit(size = 160000)
## [1] Inf

#display.brewer.pal(n = 11, name ="RdYlGn")
colorvec <- brewer.pal(n = 11, name ="RdYlGn")

col.martonne <- c("Desert" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Temperate" = colorvec[8], "Humid"= colorvec[10],"Forest"= colorvec[11])
#col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "NA" = "#D5D8DC")

col.cat <- c("Hyperarid" = "#D73027", "Arid"= "#FDAE61","Semi-arid"= "#FEE08B", "Dry subhumid" = "#ABD9E9", "Humid"= "#4575B4", "Cold" = "#DAE2F7")

1 Data

Liens utiles:

• https://www.fao.org/3/x0490E/x0490e0k.htm#annex%203.%20background%20on%20physical%20parameters%20used%20in%20evapotranspiration%20computatio

# era5vA <- read.table("era5_AI.txt") %>%
#   mutate(tas = tas-273.15, pr = mpr * 60 * 60 * 24 * 365, rsds = ssrd / (60*60*24) * 1e-6 * 60*60*24) %>% # conversion from J/m2 to J/sm2 to MJ/m2/day  
#   select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn","rsds","ea","ET0","AI","cat.AI"))
# 
# 
# write.table(era5vA, "era5vA.txt")

era5vA <- read.table("ERA5/era5_AI.txt") %>% 
  mutate(z = z, tas = t2m - 273.15, sfcWind = wind, model = "historical", source = "ERA5") %>%  #pr, ET0 already in mm/year
  select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn","ea","z","ET0","AI","cat.AI")) %>%
  filter(lm == 1)

write.table(era5vA, "era5vA.txt")

wdcvA <- read.table("wdc_AI.txt") %>% mutate(z = z, ea = vapr, ET0 = ET0) %>% # srad in MJ/m2/day, pr and ETO already in mm/year
 select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn", "ea","z","ET0","AI","cat.AI")) %>%
  filter(lm == 1)


write.table(wdcvA, "wdcvA.txt")


cmip6 <- read.table("cmip6_AI.txt")

cmip6vA <- cmip6 %>% subset(period == "1970_2000") %>%
   group_by(lon,lat,Continent, Type, Name, Acronym, lm, model, z) %>% 
   dplyr::summarise(tas = mean(tas, na.rm = T) - 273.15, pr = mean(spr, na.rm = T), # pr in mm/year
                   sfcWind = mean(sfcWind, na.rm = T), Rn = mean(Rn, na.rm = T),
                   ea = mean(ea, na.rm = T), ET0 = mean(sET0, na.rm = T), AI = mean(AI, na.rm = T)) %>% # ET0 in mm/year
  ungroup() %>%  
  mutate(source = "CMIP6")%>%
   select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn", "ea","z","ET0","AI"))


breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
#col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")
col.cat <- c("Hyperarid" = "#D73027", "Arid"= "#FDAE61","Semi-arid"= "#FEE08B", "Dry subhumid" = "#ABD9E9", "Humid"= "#4575B4", "Cold" = "#DAE2F7")

cmip6vA$cat.AI <- cmip6vA$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6vA$cat.AI[which(cmip6vA$ET0 < 400)] <- "Cold"

write.table(cmip6vA, "cmip6vA.txt")

era5vA <- read.table("figures_article/era5vA.txt")
wdcvA <- read.table("figures_article/wdcvA.txt")
cmip6vA <- read.table("figures_article/cmip6vA.txt")

2 Regions and Continents

colors_continents <- c("#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd",
                  "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf",
                  "#aec7e8", "#ffbb78", "#98df8a", "#ff9896")

ggplot(era5vA)+geom_raster(aes(x= lon, y = lat, fill = Continent))+
  scale_fill_discrete(type = colors_continents)+
  theme_void()

“Continents” to remove are: * Arctic (too little gridcells) * Atlantic (ocean) * Indian (ocean) * Pacific (ocean) * Southern (ocean)

land_mask <- raster("Masks/land_sea_mask_1degree.nc4") 
plot(land_mask)

land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land


elev <- raster("Worldclim/wc2.1_10m/wc2.1_10m_elev.tif") 
elev.df <- elev %>% projectRaster(to = land_mask) %>% as.data.frame(xy = T) %>% setNames(c("lon","lat", "z"))


ipcc_regions <- shapefile("Masks/IPCC-WGI-reference-regions-v4.shp") %>% spTransform(crs("EPSG:4326")) 

ipcc_regions.raster <- ipcc_regions %>% rasterize(land_mask)
ipcc_regions.df <- as.data.frame(ipcc_regions.raster, xy = T) %>% setNames(c("lon","lat","Continent","Type","Name","Acronym"))

ipcc_sf <- st_as_sf(ipcc_regions)

ipcc_sf$Name2 <- ipcc_sf$Name %>% str_wrap(width = 10)

ggplot(data = filter(ipcc_sf, Type %in% c("Land","Land-Ocean")))+
  geom_sf(aes(fill=Continent))+
  scale_fill_manual(values = colorvec)+
 geom_sf_text(aes(label = Name), size = 3, position = position_jitter(height = 3, width = 3, seed = 1))+
  #geom_sf_text(aes(label = Acronym), size = 2)+
  theme_void()+
  theme(legend.position = "bottom")

#ggsave("ipcc_regions.png", width = 12, height = 8, units = "cm", dpi = "retina")

3 Aridity Index

3.1 Comparison of WDC, ERA5 and CMIP6 for the reference period 1970-2000

calib <- rbind(cmip6vA, era5vA, wdcvA) %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC"))

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
#col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")
col.cat <- c("Hyperarid" = "#D73027", "Arid"= "#FDAE61","Semi-arid"= "#FEE08B", "Dry subhumid" = "#ABD9E9", "Humid"= "#4575B4", "Cold" = "#DAE2F7")

calib$cat.AI <- calib$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
calib$cat.AI[which(calib$ET0 < 400)] <- "Cold"

calib$cat.AI <- factor(calib$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

ggarrange(
  
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "Wordclim, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "ERA5, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "CMIP6, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = cat.AI))+
  geom_raster(aes(x = lon, y = lat), fill = "white")+
    scale_fill_manual(values = col.cat)+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "bottom"),

nrows = 4, ncol = 1

)

wdcfordiff <- calib %>% subset(source == "Worldclim") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.wdc")) 
era5fordiff <- calib %>% subset(source == "ERA5") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.era5")) 
cmip6fordiff <- calib %>% subset(source == "CMIP6") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cmip6")) 

diffcatAI <- wdcfordiff %>% merge(era5fordiff, by = c("lon","lat")) %>% merge(cmip6fordiff, by = c("lon","lat"))


ggplot(diffcatAI) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI.era5))+
    scale_fill_manual(values = col.cat)+
    borders()+
    labs(title = " Aridity categories ERA5, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggplot(filter(diffcatAI, cat.AI.wdc != cat.AI.era5)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI.wdc))+
    scale_fill_manual(values = col.cat)+
    borders()+
    labs(title = "Difference Worldclim/ERA5, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggplot(filter(diffcatAI, cat.AI.cmip6 != cat.AI.era5)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI.cmip6))+
    scale_fill_manual(values = col.cat)+
    borders()+
    labs(title = "Difference CMIP6/ERA5, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggplot(filter(diffcatAI, cat.AI.cmip6 != cat.AI.wdc)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI.cmip6))+
    scale_fill_manual(values = col.cat)+
    borders()+
    labs(title = "Difference CMIP6/Worldclim, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

3.2 Comparison of WDC, ERA5 and CMIP6 for the reference period 1970-2000 with country borders

calib <- rbind(cmip6vA, era5vA, wdcvA)

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
#col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")
col.cat <- c("Hyperarid" = "#D73027", "Arid"= "#FDAE61","Semi-arid"= "#FEE08B", "Dry subhumid" = "#ABD9E9", "Humid"= "#4575B4", "Cold" = "#DAE2F7")

calib$cat.AI <- calib$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
calib$cat.AI[which(calib$ET0 < 400)] <- "Cold"

calib$cat.AI <- factor(calib$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

ggarrange(
  
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "Wordclim, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "ERA5, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "CMIP6, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = cat.AI))+
  geom_raster(aes(x = lon, y = lat), fill = "white")+
    scale_fill_manual(values = col.cat)+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "top"),

nrows = 4, ncol = 1

)

ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggsave("worldclim_AIref.png", width = 12, height = 8, dpi = "retina")

ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggsave("era5_AIref.png", width = 12, height = 8, dpi = "retina")

ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggsave("cmip6_AIref.png", width = 12, height = 8, dpi = "retina")

ggarrange(
  
ggplot(data = subset(calib, source == "Worldclim" & AI < 0.05 & AI > 0)) + geom_raster(aes(x=lon, y = lat,  fill = AI))+
    scale_fill_gradient(low = "#D73027", high = "#4575B4")+
    borders(colour= "grey60")+
    labs(title = "Wordclim, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "ERA5" & AI < 0.05 & AI > 0)) + geom_raster(aes(x=lon, y = lat,  fill = AI))+
    scale_fill_gradient(low = "#D73027", high = "#4575B4")+
    borders(colour= "grey60")+
    labs(title = "ERA5, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "CMIP6" & AI < 0.05 & AI > 0)) + geom_raster(aes(x=lon, y = lat,  fill = AI))+
    scale_fill_gradient(low = "#D73027", high = "#4575B4")+
    borders(colour= "grey60")+
    labs(title = "CMIP6, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "Worldclim" & AI < 0.05 & AI > 0)) + geom_raster(aes(x = lon, y = lat, fill = AI))+
  geom_raster(aes(x = lon, y = lat), fill = "white")+
    scale_fill_gradient(low = "#D73027", high = "#4575B4")+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "top"),

nrows = 4, ncol = 1

)

3.2.1 Compare ETO and pr in WDC, ERA5, CMIP6

bpr <- c(-400, 0, 400, 800, 1200, 1600, 5000)
colvecrdbu <-  brewer.pal(n = 11, name ="RdYlBu")

col.pr <- c("(-400,0]" = colvecrdbu[5], "(0,400]"= colvecrdbu[6],"(400,800]"= colvecrdbu[7], "(800,1.2e+03]" = colvecrdbu[8], "(1.2e+03,1.6e+03]"= colvecrdbu[9], "(1.6e+03,5e+03]" = colvecrdbu[10])

calib$pr.cat <- calib$pr %>% cut(breaks = bpr) 

g.pr.wdc <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = pr.cat))+
    scale_fill_manual(values= col.pr, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual precipitation in Wordclim", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

g.pr.era5 <- ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = pr.cat))+
      scale_fill_manual(values= col.pr, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual precipitation in ERA5", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

g.pr.cmip6 <- ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = pr.cat))+
    scale_fill_manual(values= col.pr, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual precipitation CMIP6", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

l.pr <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = pr.cat))+
    scale_fill_manual(values= col.pr, na.translate = F)+
  geom_raster(aes(x = lon, y = lat), fill = "white")+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "top", legend.key.width = unit(2, "cm"))

#)

bpr <- c(-400, 0, 400, 800, 1200, 1600, 5000)
colvecrdbu <-  brewer.pal(n = 11, name ="RdYlBu")

col.et0 <- c("(-400,0]" = colvecrdbu[8], "(0,400]"= colvecrdbu[7],"(400,800]"= colvecrdbu[6], "(800,1.2e+03]" = colvecrdbu[5], "(1.2e+03,1.6e+03]"= colvecrdbu[4], "(1.6e+03,5e+03]" = colvecrdbu[3])

calib$et0.cat <- calib$ET0 %>% cut(breaks = bpr) 

g.et0.wdc <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = et0.cat))+
   scale_fill_manual(values= col.et0, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual evapotranspiration Wordclim", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

g.et0.era5 <- ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = et0.cat))+
   scale_fill_manual(values= col.et0, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual evapotranspiration ERA5", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

g.et0.cmip6 <- ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = et0.cat))+
   scale_fill_manual(values= col.et0, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual evapotranspiration CMIP6", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

l.et0 <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = et0.cat))+
   scale_fill_manual(values= col.et0, na.translate = F)+
    geom_raster(aes(x = lon, y = lat), fill = "white")+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "top", legend.key.width = unit(2, "cm"))

ggarrange(
  ggarrange(plotlist = list(g.pr.wdc,g.pr.era5, g.pr.cmip6, l.pr), nrow = 4, ncol = 1),
  ggarrange(plotlist = list(g.et0.wdc,g.et0.era5, g.et0.cmip6, l.et0), nrow = 4, ncol = 1),
  ncol = 2)

3.3 Comparison of WDC, ERA5 and CMIP6 for the reference period 1970-2000 with Antarctica

3.3.1 Map aridity categories for the 3 databases

calib <- rbind(cmip6vA, era5vA, wdcvA)

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
# col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")
col.cat <- c("Hyperarid" = "#D73027", "Arid"= "#FDAE61","Semi-arid"= "#FEE08B", "Dry subhumid" = "#ABD9E9", "Humid"= "#4575B4", "Cold" = "#DAE2F7")

calib$cat.AI <- calib$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
calib$cat.AI[which(calib$ET0 < 400)] <- "Cold"
calib$cat.AI <- factor(calib$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

ggarrange(
  
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "Wordclim, 1970-2000", fill = "")+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "ERA5, 1970-2000", fill = "")+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "CMIP6, 1970-2000", fill = "")+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = cat.AI))+
  geom_raster(aes(x = lon, y = lat), fill = "white")+
    scale_fill_manual(values = col.cat)+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "bottom"),

nrows = 4, ncol = 1

)

3.3.2 Proportion aridity category by subregion, continent et database

calib$same.cat <- NA
calib$cat.ref <- NA

for(i in 1:nrow(calib)){
  loni = calib$lon[i]
  lati = calib$lat[i]
  catrefi <- subset(calib, lon == loni & lat == lati & source == "CMIP6")$cat.AI %>% as.character()
  cati <- calib$cat.AI[i] %>% as.character()
  calib$same.cat[i] <- ifelse(cati == catrefi, "y",
                              ifelse(cati!= catrefi, "n", 
                                     NA))
  calib$cat.ref[i] <- catrefi
}

write.table(calib, "figures_article/calib.txt")

# table_calib_c <- calib %>% select(c("Continent","source","cat.AI")) %>%
#   reshape2::melt(id.vars = c("Continent","source")) %>%
#   group_by(Continent,source,value) %>%
#   summarise(count = n()) %>% mutate(perc = round(count/sum(count)*100,1)) %>%
#   ungroup()
# 
# table_calib_c$value[is.na(table_calib_c$value) == T] <- "NA"
# 
# ggplot(table_calib_c)+
#   geom_col(aes(x = Continent, y = perc, fill = value, position = "fill"))+
#   scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","NA"))+
#   theme_minimal()+
#   theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))+
#   facet_grid(rows = vars(source))+
#   labs(x = "Continent", y = "%", fill = "Aridity\ncategory", y = "Proportion of aridity category")
# 


table_calib <- calib %>% select(c("Continent","Name","Acronym","source","cat.AI")) %>%
  reshape2::melt(id.vars = c("Continent","Name","Acronym","source")) %>%
  filter(!is.na(value)) %>%
  group_by(Continent,Name,Acronym,source,value) %>%
  summarise(count = n()) %>% mutate(perc = round(count/sum(count)*100,1)) %>%
  ungroup() %>% filter(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC"))

table_calib$value[is.na(table_calib$value) == T] <- "NA"
table_calib$Continent[which(table_calib$Continent == "EUROPE-AFRICA")] <- "EUROPE"
table_calib$Continent[which(table_calib$Continent == "CENTRAL-AMERICA")] <- "CENTRAL\nAMERICA"
table_calib$Continent[which(table_calib$Continent == "NORTH-AMERICA")] <- "NORTH\nAMERICA"
table_calib$Continent[which(table_calib$Continent == "SOUTH-AMERICA")] <- "SOUTH\nAMERICA"

table_calib$Continent <- factor(table_calib$Continent, levels = c("AFRICA","EUROPE","ASIA","OCEANIA","NORTH\nAMERICA","CENTRAL\nAMERICA","SOUTH\nAMERICA","POLAR","NA"))

table_calib$value <- factor(table_calib$value, levels = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","Cold"))

ggplot(table_calib)+
  geom_col(aes(x = Acronym, y = perc, fill = value, position = "fill"))+
  scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","Cold"))+
  theme_minimal(base_size = 12)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), 
        legend.position = "bottom", 
        strip.text.x = element_text(size = 6))+
  facet_grid(source ~ Continent, switch = "y", scales = "free_x", space = "free")+
  labs(x = "", y = "%", fill = "Aridity\ncategory")

ggplot(table_calib)+
  geom_col(aes(x = Acronym, y = count, fill = value, position = "fill"))+
  scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","Cold"))+
  scale_y_continuous(position = "right")+
  theme_minimal(base_size = 12)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), 
        legend.position = "bottom", 
        strip.text.x = element_text(size = 6))+
  facet_grid(source ~ Continent, switch = "y", scales = "free_x", space = "free")+
  labs(x = "", y = "Count", fill = "Aridity\ncategory")

3.3.3 Pie chart

library(ggrepel)

pie_calib <- table_calib %>% group_by(source, value) %>% summarise(gridcells = sum(count)) %>% mutate(percent = round(gridcells/sum(gridcells)*100, 1), lab.y = (rev(cumsum(rev(percent)))) - percent*0.5)


ggplot(pie_calib)+geom_bar(aes(x = "", y = percent, fill = value), width = 1, stat = "identity")+
  coord_polar("y", start = 0)+
  scale_fill_manual(values = col.cat, na.translate = F)+
  facet_grid(cols = vars(source))+
  geom_label_repel(aes(x = "", y = lab.y, label = percent), nudge_x = 0.5)+
  labs(fill = "")+
  theme_void()+theme(legend.position = "bottom")

pie_calib <- rbind(cmip6vA, era5vA, wdcvA)

pie_calib$cat.AI <- factor(pie_calib$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold")) 

pie_calib <- pie_calib %>% 
  subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  select(c("source","cat.AI")) %>%
   table() %>% as.data.frame() %>%
   group_by(source) %>%
    mutate(perc = round(Freq/sum(Freq)*100, 1),
         lab.y = rev(cumsum(rev(Freq))) - Freq*0.5) 

ggplot(pie_calib)+geom_bar(aes(x = "", y = Freq, fill = cat.AI), width = 1, stat = "identity")+
  coord_polar("y", start = 0)+
  scale_fill_manual(values = col.cat, na.translate = F)+
  facet_grid(cols = vars(source))+
  geom_label_repel(aes(x = "", y = lab.y, label = perc), nudge_x = 0.5)+
  labs(fill = "")+
  theme_void()+theme(legend.position = "bottom")

# cmip6 <- read.table("cmip6_AI.txt")
# 
# cmip6.ref <- cmip6 %>% filter(period == "1970_2000") %>% select(c("lon","lat","model","period","source","AI","spr","t","sET0")) %>% setNames(c("lon","lat","model","period","source","AI.ref","spr.ref","t.ref","ET0.ref"))
# cmip6.ref$cat.AI <- cmip6.ref$AI.ref %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
# cmip6.ref$cat.AI[which(cmip6.ref$ET0.ref < 400)] <- "Cold"
# 
# write.table(cmip6.ref, "figures_article/cmip6.ref.txt")

cmip6.ref <- read.table("figures_article/cmip6.ref.txt")

wdc <- filter(calib, source == "Worldclim") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.wdc"))
era5 <- filter(calib, source == "ERA5") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.era5"))
cmip6.cal <- filter(calib, source == "CMIP6") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cmip6"))

casesm2 <- filter(cmip6.ref, source == "CAS-ESM2") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.casesm2"))
cesm <- filter(cmip6.ref, source == "CESM") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cesm"))
cmcc <- filter(cmip6.ref, source == "CMCC") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cmcc"))
cmccesm2 <- filter(cmip6.ref, source == "CMCC-ESM2") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cmccesm2"))
cnrm <- filter(cmip6.ref, source == "CNRM") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cnrm"))
ecearth3 <- filter(cmip6.ref, source == "EC-Earth3") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.ecearth3"))
fgoals <- filter(cmip6.ref, source == "FGOALS") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.fgoals"))
gfdlesm4 <- filter(cmip6.ref, source == "GFDL-ESM4") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.gfdlesm4"))
inm <- filter(cmip6.ref, source == "INM") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.inm"))
inmcm5 <- filter(cmip6.ref, source == "INM-CM5") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.inmcm5"))
mpi <- filter(cmip6.ref, source == "MPI") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.mpi"))
mri <- filter(cmip6.ref, source == "MRI") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.mri"))
noresm <- filter(cmip6.ref, source == "NorESM-2-MM") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.noresm"))

all <- merge(wdc, era5, by = c("lon","lat")) %>% 
  merge(cmip6.cal, by = c("lon","lat")) %>% 
  merge(casesm2, by = c("lon","lat")) %>%
  merge(cesm, by = c("lon","lat")) %>%
  merge(cmcc, by = c("lon","lat")) %>%
  merge(cmccesm2, by = c("lon","lat")) %>%
  merge(cnrm, by = c("lon","lat")) %>%
  merge(ecearth3, by = c("lon","lat")) %>%
  merge(fgoals, by = c("lon","lat")) %>%
  merge(gfdlesm4, by = c("lon","lat")) %>%
  merge(inm, by = c("lon","lat")) %>%
  merge(inmcm5, by = c("lon","lat")) %>%
  merge(mpi, by = c("lon","lat")) %>%
  merge(mri, by = c("lon","lat")) %>%
  merge(noresm, by = c("lon","lat"))

all$cmip6.wdc <- with(all, ifelse(cat.AI.cmip6 == cat.AI.wdc, "same","diff"))
all$cmip6.era5 <- with(all, ifelse(cat.AI.cmip6 == cat.AI.era5, "same","diff"))
all$wdc.era5 <- with(all, ifelse(cat.AI.wdc == cat.AI.era5, "same","diff"))


all$casesm2.wdc <- with(all, ifelse(cat.AI.casesm2 == cat.AI.wdc, "same","diff"))
all$cesm.wdc <- with(all, ifelse(cat.AI.cesm == cat.AI.wdc, "same","diff"))
all$cmcc.wdc <- with(all, ifelse(cat.AI.cmcc == cat.AI.wdc, "same","diff"))
all$cmccesm2.wdc <- with(all, ifelse(cat.AI.cmccesm2 == cat.AI.wdc, "same","diff"))
all$cnrm.wdc <- with(all, ifelse(cat.AI.cnrm == cat.AI.wdc, "same","diff"))
all$ecearth3.wdc <- with(all, ifelse(cat.AI.ecearth3 == cat.AI.wdc, "same","diff"))
all$fgoals.wdc <- with(all, ifelse(cat.AI.fgoals == cat.AI.wdc, "same","diff"))
all$gfdlesm4.wdc <- with(all, ifelse(cat.AI.gfdlesm4 == cat.AI.wdc, "same","diff"))
all$inm.wdc <- with(all, ifelse(cat.AI.inm == cat.AI.wdc, "same","diff"))
all$inmcm5.wdc <- with(all, ifelse(cat.AI.inmcm5 == cat.AI.wdc, "same","diff"))
all$mpi.wdc <- with(all, ifelse(cat.AI.mpi == cat.AI.wdc, "same","diff"))
all$mri.wdc <- with(all, ifelse(cat.AI.mri == cat.AI.wdc, "same","diff"))
all$noresm.wdc <- with(all, ifelse(cat.AI.noresm == cat.AI.wdc, "same","diff"))

df1 <- table(all$cmip6.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CMIP6-WDC"))
df2 <- table(all$casesm2.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CASESM2-WDC"))
df3 <- table(all$cesm.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CESM-WDC"))
df4 <- table(all$cmcc.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CMCC-WDC"))
df5 <- table(all$cmccesm2.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CMCCESM2-WDC"))
df6 <- table(all$cnrm.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CNRM-WDC"))
df7 <- table(all$ecearth3.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "ECEARTH3-WDC"))
df8 <- table(all$fgoals.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "FGOALS-WDC"))
df9 <- table(all$gfdlesm4.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "GFDLESM4-WDC"))
df10 <- table(all$inm.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "INM-WDC"))
df11 <- table(all$inmcm5.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "INMCM5-WDC"))
df12 <- table(all$mpi.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "MPI-WDC"))
df13 <- table(all$mri.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "MRI-WDC"))
df14 <- table(all$noresm.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "NORESM-WDC"))

df0 <- merge(df1,df2, by = "diffornot") %>% merge(df3, by = "diffornot") %>% 
  merge(df4, by = "diffornot") %>% merge(df5, by = "diffornot") %>% merge(df6, by = "diffornot") %>% 
  merge(df7, by = "diffornot") %>% merge(df8, by = "diffornot") %>% merge(df9, by = "diffornot") %>%
  merge(df10, by = "diffornot") %>% merge(df11, by = "diffornot") %>% merge(df12, by = "diffornot") %>%
  merge(df13, by = "diffornot") %>% merge(df14, by = "diffornot")

prop.diff.wdc <- df0[1,]/(df0[1,] + df0[2,])*100      
prop.diff.wdc[1,1] <- "Percent diff"
                                                      
df.prop.wdc <- rbind(df0, prop.diff.wdc)

kable(df.prop.wdc, caption = "Percent of gridcells in which the aridity category differs from Worldclim", digits = 1) %>%
  kable_styling(bootstrap_options = "bordered")

Percent of gridcells in which the aridity category differs from Worldclim

diffornot	CMIP6-WDC	CASESM2-WDC	CESM-WDC	CMCC-WDC	CMCCESM2-WDC	CNRM-WDC	ECEARTH3-WDC	FGOALS-WDC	GFDLESM4-WDC	INM-WDC	INMCM5-WDC	MPI-WDC	MRI-WDC	NORESM-WDC
diff	3273.0	4275.0	3557.0	4006.0	6503	4594.0	3691	5083.0	3531.0	4323.0	4318.0	4276.0	3349.0	3843.0
same	18418.0	17407.0	18125.0	17676.0	15179	17097.0	18000	16608.0	18160.0	17359.0	17364.0	17415.0	18333.0	17839.0
Percent diff	15.1	19.7	16.4	18.5	30	21.2	17	23.4	16.3	19.9	19.9	19.7	15.4	17.7

all$casesm2.era5 <- with(all, ifelse(cat.AI.casesm2 == cat.AI.era5, "same","diff"))
all$cesm.era5 <- with(all, ifelse(cat.AI.cesm == cat.AI.era5, "same","diff"))
all$cmcc.era5 <- with(all, ifelse(cat.AI.cmcc == cat.AI.era5, "same","diff"))
all$cmccesm2.era5 <- with(all, ifelse(cat.AI.cmccesm2 == cat.AI.era5, "same","diff"))
all$cnrm.era5 <- with(all, ifelse(cat.AI.cnrm == cat.AI.era5, "same","diff"))
all$ecearth3.era5 <- with(all, ifelse(cat.AI.ecearth3 == cat.AI.era5, "same","diff"))
all$fgoals.era5 <- with(all, ifelse(cat.AI.fgoals == cat.AI.era5, "same","diff"))
all$gfdlesm4.era5 <- with(all, ifelse(cat.AI.gfdlesm4 == cat.AI.era5, "same","diff"))
all$inm.era5 <- with(all, ifelse(cat.AI.inm == cat.AI.era5, "same","diff"))
all$inmcm5.era5 <- with(all, ifelse(cat.AI.inmcm5 == cat.AI.era5, "same","diff"))
all$mpi.era5 <- with(all, ifelse(cat.AI.mpi == cat.AI.era5, "same","diff"))
all$mri.era5 <- with(all, ifelse(cat.AI.mri == cat.AI.era5, "same","diff"))
all$noresm.era5 <- with(all, ifelse(cat.AI.noresm == cat.AI.era5, "same","diff"))

df1 <- table(all$cmip6.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CMIP6-ERA5"))
df2 <- table(all$casesm2.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CASESM2-ERA5"))
df3 <- table(all$cesm.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CESM-ERA5"))
df4 <- table(all$cmcc.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CMCC-ERA5"))
df5 <- table(all$cmccesm2.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CMCCESM2-ERA5"))
df6 <- table(all$cnrm.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CNRM-ERA5"))
df7 <- table(all$ecearth3.era5) %>% as.data.frame() %>% setNames(c("diffornot", "ECEARTH3-ERA5"))
df8 <- table(all$fgoals.era5) %>% as.data.frame() %>% setNames(c("diffornot", "FGOALS-ERA5"))
df9 <- table(all$gfdlesm4.era5) %>% as.data.frame() %>% setNames(c("diffornot", "GFDLESM4-ERA5"))
df10 <- table(all$inm.era5) %>% as.data.frame() %>% setNames(c("diffornot", "INM-ERA5"))
df11 <- table(all$inmcm5.era5) %>% as.data.frame() %>% setNames(c("diffornot", "INMCM5-ERA5"))
df12 <- table(all$mpi.era5) %>% as.data.frame() %>% setNames(c("diffornot", "MPI-ERA5"))
df13 <- table(all$mri.era5) %>% as.data.frame() %>% setNames(c("diffornot", "MRI-ERA5"))
df14 <- table(all$noresm.era5) %>% as.data.frame() %>% setNames(c("diffornot", "NORESM-ERA5"))

df0 <- merge(df1,df2, by = "diffornot") %>% merge(df3, by = "diffornot") %>% 
  merge(df4, by = "diffornot") %>% merge(df5, by = "diffornot") %>% merge(df6, by = "diffornot") %>% 
  merge(df7, by = "diffornot") %>% merge(df8, by = "diffornot") %>% merge(df9, by = "diffornot") %>%
  merge(df10, by = "diffornot") %>% merge(df11, by = "diffornot") %>% merge(df12, by = "diffornot") %>%
  merge(df13, by = "diffornot") %>% merge(df14, by = "diffornot")

prop.diff.era5 <- df0[1,]/(df0[1,] + df0[2,])*100      
prop.diff.era5[1,1] <- "Percent diff"
                                                      
df.prop.era5 <- rbind(df0, prop.diff.era5)   

kable(df.prop.era5, caption = "Percent of gridcells in which the aridity category differs from Worldclim", digits = 1) %>%
  kable_styling(bootstrap_options = "bordered")

Percent of gridcells in which the aridity category differs from Worldclim

diffornot	CMIP6-ERA5	CASESM2-ERA5	CESM-ERA5	CMCC-ERA5	CMCCESM2-ERA5	CNRM-ERA5	ECEARTH3-ERA5	FGOALS-ERA5	GFDLESM4-ERA5	INM-ERA5	INMCM5-ERA5	MPI-ERA5	MRI-ERA5	NORESM-ERA5
diff	3197.0	3934.0	3451.0	3819.0	6969.0	3859.0	3629.0	5402.0	3860.0	4189.0	4015.0	3633.0	2910.0	3683
same	18494.0	17748.0	18231.0	17863.0	14713.0	17832.0	18062.0	16289.0	17831.0	17493.0	17667.0	18058.0	18772.0	17999
Percent diff	14.7	18.1	15.9	17.6	32.1	17.8	16.7	24.9	17.8	19.3	18.5	16.7	13.4	17

df <- table(all$wdc.era5) %>% as.data.frame() %>% setNames(c("diffornot","ERA5-WDC"))
percent.diff <- df[1,] / (df[1,]+df[2,]) * 100
percent.diff[1,1] <- "Percent diff"
df.prop.wdc.era5 <- rbind(df, percent.diff)

kable(df.prop.wdc.era5, caption = "Percent of gridcells in which the aridity category differs between Worldclim and ERA5", digits = 1) %>% 
  kable_styling(bootstrap_options = "bordered")

Percent of gridcells in which the aridity category differs between Worldclim and ERA5

diffornot	ERA5-WDC
diff	2850.0
same	18841.0
Percent diff	13.1

4 Number of inhabitants by region

https://earthdata.nasa.gov/data/catalog/sedac-ciesin-sedac-gpwv4-popcount-r11-4.11

density <- raster("Population_density/gpw_v4_population_density_rev11_2020_1_deg.tif")
plot(density)

popcount <- raster("Population_density/gpw_v4_population_count_rev11_2020_1_deg.tif")
plot(popcount)

land_mask <- raster("Masks/land_sea_mask_1degree.nc4") 
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land

popcount.df <- popcount %>%
  projectRaster(to = land_mask) %>% 
  as.data.frame(xy = T) %>% setNames(c("lon","lat","pop")) %>% 
  merge(land_mask.df, by = c("lon","lat")) %>% filter(lm == 1)

era5vApop <- merge(era5vA, popcount.df, by = c("lon","lat")) 

tab.pop.by.cat <- era5vApop %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(cat.AI, Continent) %>%
  summarise(pop = sum(pop)) 

write.csv(tab.pop.by.cat, "population.count.era5.csv")

knitr::kable(tab.pop.by.cat)

cat.AI	Continent	pop
Arid	AFRICA	9.407395e+07
Arid	ASIA	2.608739e+08
Arid	CENTRAL-AMERICA	1.102838e+07
Arid	EUROPE	3.523270e+05
Arid	EUROPE-AFRICA	3.354203e+07
Arid	NORTH-AMERICA	8.275538e+05
Arid	OCEANIA	2.483547e+05
Arid	SOUTH-AMERICA	2.714742e+06
Cold	ASIA	9.005723e+06
Cold	EUROPE	6.822776e+06
Cold	NORTH-AMERICA	8.555420e+05
Cold	POLAR	NA
Cold	SOUTH-AMERICA	6.222267e+04
Dry subhumid	AFRICA	1.131375e+08
Dry subhumid	ASIA	3.724373e+08
Dry subhumid	CENTRAL-AMERICA	2.847546e+07
Dry subhumid	EUROPE	5.195230e+06
Dry subhumid	EUROPE-AFRICA	5.933504e+07
Dry subhumid	NORTH-AMERICA	9.392606e+06
Dry subhumid	OCEANIA	7.027896e+05
Dry subhumid	SOUTH-AMERICA	1.297947e+07
Humid	AFRICA	6.854010e+08
Humid	ASIA	3.055659e+09
Humid	CENTRAL-AMERICA	1.462475e+08
Humid	EUROPE	4.999673e+08
Humid	EUROPE-AFRICA	1.452014e+08
Humid	NORTH-AMERICA	2.915502e+08
Humid	OCEANIA	1.529713e+07
Humid	POLAR	4.329376e+02
Humid	SOUTH-AMERICA	3.724919e+08
Hyperarid	AFRICA	5.328568e+07
Hyperarid	ASIA	1.970440e+07
Hyperarid	EUROPE-AFRICA	4.146658e+07
Hyperarid	SOUTH-AMERICA	5.088920e+05
Semi-arid	AFRICA	2.154724e+08
Semi-arid	ASIA	5.713649e+08
Semi-arid	CENTRAL-AMERICA	3.973955e+07
Semi-arid	EUROPE	1.841667e+06
Semi-arid	EUROPE-AFRICA	8.619193e+07
Semi-arid	NORTH-AMERICA	2.806716e+07
Semi-arid	OCEANIA	3.403064e+06
Semi-arid	SOUTH-AMERICA	2.872740e+07
NA	AFRICA	3.321493e+04
NA	ASIA	1.762974e+07
NA	CENTRAL-AMERICA	7.132541e+05
NA	EUROPE	1.481967e+05
NA	EUROPE-AFRICA	1.876606e+06
NA	NORTH-AMERICA	1.403549e+03
NA	OCEANIA	3.896168e+04
NA	POLAR	NA
NA	SOUTH-AMERICA	2.214226e+05

5 Comparison of datasets

vars <- c("lon","lat","Continent","Name","Acronym","tas","pr","sfcWind","Rn","ea","ET0","AI","cat.AI")
merge.vars <- c("lon","lat","Continent","Name","Acronym")
calib.wide <- select(era5vA, vars) %>%
               merge(select(wdcvA, vars), 
                     by = merge.vars) %>%
          merge(select(cmip6vA, vars), 
                by = merge.vars) %>%
          setNames(c(merge.vars,
                     "tas.era5","pr.era5", "sfcWind.era5","Rn.era5","ea.era5","ET0.era5","AI.era5","cat.AI.era5",
                     "tas.wdc","pr.wdc", "sfcWind.wdc","Rn.wdc","ea.wdc","ET0.wdc","AI.wdc","cat.AI.wdc",
                     "tas.cmip6","pr.cmip6", "sfcWind.cmip6","Rn.cmip6","ea.cmip6","ET0.cmip6","AI.cmip6","cat.AI.cmip6"))

5.1 Scatterplots with Antarctica

5.1.1 ERA5 compared to WDC

Rn is computed differently for ERA5 and WDC, hence the discrepancy for low values.

5.1.2 ERA5 and WDC

ggarrange(

ggplot(calib.wide,aes(x = tas.wdc))+geom_point(aes( y = tas.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, ERA5", col = "Dataset")+
  scale_color_manual(values = "#A50026")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = "#A50026")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = "#A50026")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = "#A50026")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = "#A50026")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = "#A50026")+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

knitr::kable(data.frame(precipitation = cor(calib.wide$pr.wdc, calib.wide$pr.era5, use = "pairwise.complete.obs"),
                        temperature = cor(calib.wide$tas.wdc, calib.wide$tas.era5, use = "pairwise.complete.obs"),
                        wind = cor(calib.wide$sfcWind.wdc, calib.wide$sfcWind.era5, use = "pairwise.complete.obs"),
                        Rn = cor(calib.wide$Rn.wdc, calib.wide$Rn.era5, use = "pairwise.complete.obs"),
                        ea = cor(calib.wide$ea.wdc, calib.wide$ea.era5, use = "pairwise.complete.obs"),
                        ET0 = cor(calib.wide$ET0.wdc, calib.wide$ET0.era5, use = "pairwise.complete.obs")), 
             digits = 2, caption = "Correlation coefficients between ERA5 and Worldclim") %>%
  kableExtra::kable_styling(bootstrap_options = "bordered")

Correlation coefficients between ERA5 and Worldclim

precipitation	temperature	wind	Rn	ea	ET0
0.89	0.99	0.68	0.67	1	0.98

ggplot(calib.wide) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none")


ggarrange(

ggplot(calib.wide) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

5.1.3 CMIP6 and WDC

ggarrange(

ggplot(calib.wide,aes(x = tas.wdc))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, CMIP6")+
  scale_color_manual(values = "#FDAE61")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.wdc))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, CMIP6")+
  scale_color_manual(values = "#FDAE61")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.wdc))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, CMIP6")+
  scale_color_manual(values = "#FDAE61")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.wdc))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, CMIP6")+
   scale_color_manual(values = "#FDAE61")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.wdc))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, CMIP6")+
   scale_color_manual(values = "#FDAE61")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.wdc))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, CMIP6")+
   scale_color_manual(values = "#FDAE61")+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2
)

5.1.4 CMIP6 and ERA5

ggarrange(
  
ggplot(calib.wide,aes(x = tas.era5))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, ERA5", y = "Temperature in °C, CMIP6", col = "Dataset")+
  scale_color_manual(values = "#4575B4")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.era5))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, ERA5", y = "Annual precipitation, mm, CMIP6", col = "Dataset")+
  scale_color_manual(values = "#4575B4")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.era5))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s, CMIP6", col = "Dataset")+
  scale_color_manual(values = "#4575B4")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.era5))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
   stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, ERA5, °C", y = "Rn - G, W/MJ/m2/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = "#4575B4")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.era5))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, ERA5", y = "ea, kPa, CMIP6", col = "Dataset")+
   scale_color_manual(values = "#4575B4")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.era5))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, ERA5", y = "ET0, mm/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = "#4575B4")+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2)

5.2 Scatterplots without Antarctica

calib.noant <- calib.wide %>% filter(lat > -55)

5.2.1 Boxplots

library(ggtext)

ggarrange(
  
ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = pr.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = pr.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = pr.cmip6, fill = "CMIP6"))+
  labs(x = "Mean annual precipitation, mm", y = "", fill = "")+
  xlim(0,6000)+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = tas.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = tas.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = tas.cmip6, fill = "CMIP6"))+
  labs(x = "Surface temperature, °C", y = "", fill = "")+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = sfcWind.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = sfcWind.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = sfcWind.cmip6, fill = "CMIP6"))+
  labs(x = expression("Surface wind speed, " * m~s^{-1}), y = "", fill = "")+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = Rn.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = Rn.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = Rn.cmip6, fill = "CMIP6"))+
  labs(x = expression("Rn - G, " * MJ~m^{-2}~day^{-1}), y = "", fill = "")+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = ea.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = ea.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = ea.cmip6, fill = "CMIP6"))+
  labs(x = "ea, kPa", y = "", fill = "")+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = ET0.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = ET0.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = ET0.cmip6, fill = "CMIP6"))+
  labs(x = expression("ET0, " * mm~day^{-1}), y = "", fill = "")+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ncol = 3, nrow = 2, common.legend = T, legend = "bottom"
)

library(SpatialPack)


coords <- calib.noant[,c("lon","lat")]

df.era5.wdc <- with(calib.noant, 
                    data.frame(precipitation = c(cor(pr.wdc, pr.era5, use = "pairwise.complete.obs"), cor.test(pr.wdc, pr.era5, use = "pairwise.complete.obs")[[3]], cor(pr.wdc, pr.era5, use = "pairwise.complete.obs")^2, modified.ttest(pr.wdc, pr.era5, coords)$corr, modified.ttest(pr.wdc, pr.era5, coords)$p.value),
                        temperature = c(cor(tas.wdc, tas.era5, use = "pairwise.complete.obs"), cor.test(tas.wdc, tas.era5, use = "pairwise.complete.obs")[[3]], cor(tas.wdc, tas.era5, use = "pairwise.complete.obs")^2, modified.ttest(tas.wdc, tas.era5, coords)$corr, modified.ttest(tas.wdc, tas.era5, coords)$p.value),
                        wind = c(cor(sfcWind.wdc, sfcWind.era5, use = "pairwise.complete.obs"), cor.test(sfcWind.wdc, sfcWind.era5, use = "pairwise.complete.obs")[[3]], cor(sfcWind.wdc, sfcWind.era5, use = "pairwise.complete.obs")^2, modified.ttest(sfcWind.wdc, sfcWind.era5, coords)$corr, modified.ttest(sfcWind.wdc, sfcWind.era5, coords)$p.value),
                        Rn = c(cor(Rn.wdc, Rn.era5, use = "pairwise.complete.obs"), cor.test(Rn.wdc, Rn.era5, use = "pairwise.complete.obs")[[3]], cor(Rn.wdc, Rn.era5, use = "pairwise.complete.obs")^2, modified.ttest(Rn.wdc, Rn.era5, coords)$corr, modified.ttest(Rn.wdc, Rn.era5, coords)$p.value),
                        ea = c(cor(ea.wdc, ea.era5, use = "pairwise.complete.obs"), cor.test(ea.wdc, ea.era5, use = "pairwise.complete.obs")[[3]], cor(ea.wdc, ea.era5, use = "pairwise.complete.obs")^2, modified.ttest(ea.wdc, ea.era5, coords)$corr, modified.ttest(ea.wdc, ea.era5, coords)$p.value),
                        ET0 = c(cor(ET0.wdc, ET0.era5, use = "pairwise.complete.obs"), cor.test(ET0.wdc, ET0.era5, use = "pairwise.complete.obs")[[3]], cor(ET0.wdc, ET0.era5, use = "pairwise.complete.obs")^2, modified.ttest(ET0.wdc, ET0.era5, coords)$corr, modified.ttest(ET0.wdc, ET0.era5, coords)$p.value)))

df.era5.cmip6 <- with(calib.noant, 
                    data.frame(precipitation = c(cor(pr.cmip6, pr.era5, use = "pairwise.complete.obs"), cor.test(pr.cmip6, pr.era5, use = "pairwise.complete.obs")[[3]], cor(pr.cmip6, pr.era5, use = "pairwise.complete.obs")^2, modified.ttest(pr.cmip6, pr.era5, coords)$corr, modified.ttest(pr.cmip6, pr.era5, coords)$p.value),
                        temperature = c(cor(tas.cmip6, tas.era5, use = "pairwise.complete.obs"), cor.test(tas.cmip6, tas.era5, use = "pairwise.complete.obs")[[3]], cor(tas.cmip6, tas.era5, use = "pairwise.complete.obs")^2, modified.ttest(tas.cmip6, tas.era5, coords)$corr, modified.ttest(tas.cmip6, tas.era5, coords)$p.value),
                        wind = c(cor(sfcWind.cmip6, sfcWind.era5, use = "pairwise.complete.obs"), cor.test(sfcWind.cmip6, sfcWind.era5, use = "pairwise.complete.obs")[[3]], cor(sfcWind.cmip6, sfcWind.era5, use = "pairwise.complete.obs")^2, modified.ttest(sfcWind.cmip6, sfcWind.era5, coords)$corr, modified.ttest(sfcWind.cmip6, sfcWind.era5, coords)$p.value),
                        Rn = c(cor(Rn.cmip6, Rn.era5, use = "pairwise.complete.obs"), cor.test(Rn.cmip6, Rn.era5, use = "pairwise.complete.obs")[[3]], cor(Rn.cmip6, Rn.era5, use = "pairwise.complete.obs")^2, modified.ttest(Rn.cmip6, Rn.era5, coords)$corr, modified.ttest(Rn.cmip6, Rn.era5, coords)$p.value),
                        ea = c(cor(ea.cmip6, ea.era5, use = "pairwise.complete.obs"), cor.test(ea.cmip6, ea.era5, use = "pairwise.complete.obs")[[3]], cor(ea.cmip6, ea.era5, use = "pairwise.complete.obs")^2, modified.ttest(ea.cmip6, ea.era5, coords)$corr, modified.ttest(ea.cmip6, ea.era5, coords)$p.value),
                        ET0 = c(cor(ET0.cmip6, ET0.era5, use = "pairwise.complete.obs"), cor.test(ET0.cmip6, ET0.era5, use = "pairwise.complete.obs")[[3]], cor(ET0.cmip6, ET0.era5, use = "pairwise.complete.obs")^2, modified.ttest(ET0.cmip6, ET0.era5, coords)$corr, modified.ttest(ET0.cmip6, ET0.era5, coords)$p.value)))

df.wdc.cmip6 <- with(calib.noant, 
                    data.frame(precipitation = c(cor(pr.cmip6, pr.wdc, use = "pairwise.complete.obs"), cor.test(pr.cmip6, pr.wdc, use = "pairwise.complete.obs")[[3]], cor(pr.cmip6, pr.wdc, use = "pairwise.complete.obs")^2, modified.ttest(pr.cmip6, pr.wdc, coords)$corr, modified.ttest(pr.cmip6, pr.wdc, coords)$p.value),
                        temperature = c(cor(tas.cmip6, tas.wdc, use = "pairwise.complete.obs"), cor.test(tas.cmip6, tas.wdc, use = "pairwise.complete.obs")[[3]], cor(tas.cmip6, tas.wdc, use = "pairwise.complete.obs")^2, modified.ttest(tas.cmip6, tas.wdc, coords)$corr, modified.ttest(tas.cmip6, tas.wdc, coords)$p.value),
                        wind = c(cor(sfcWind.cmip6, sfcWind.wdc, use = "pairwise.complete.obs"), cor.test(sfcWind.cmip6, sfcWind.wdc, use = "pairwise.complete.obs")[[3]], cor(sfcWind.cmip6, sfcWind.wdc, use = "pairwise.complete.obs")^2, modified.ttest(sfcWind.cmip6, sfcWind.wdc, coords)$corr, modified.ttest(sfcWind.cmip6, sfcWind.wdc, coords)$p.value),
                        Rn = c(cor(Rn.cmip6, Rn.wdc, use = "pairwise.complete.obs"), cor.test(Rn.cmip6, Rn.wdc, use = "pairwise.complete.obs")[[3]], cor(Rn.cmip6, Rn.wdc, use = "pairwise.complete.obs")^2, modified.ttest(Rn.cmip6, Rn.wdc, coords)$corr, modified.ttest(Rn.cmip6, Rn.wdc, coords)$p.value),
                        ea = c(cor(ea.cmip6, ea.wdc, use = "pairwise.complete.obs"), cor.test(ea.cmip6, ea.wdc, use = "pairwise.complete.obs")[[3]], cor(ea.cmip6, ea.wdc, use = "pairwise.complete.obs")^2, modified.ttest(ea.cmip6, ea.wdc, coords)$corr, modified.ttest(ea.cmip6, ea.wdc, coords)$p.value),
                        ET0 = c(cor(ET0.cmip6, ET0.wdc, use = "pairwise.complete.obs"), cor.test(ET0.cmip6, ET0.wdc, use = "pairwise.complete.obs")[[3]], cor(ET0.cmip6, ET0.wdc, use = "pairwise.complete.obs")^2, modified.ttest(ET0.cmip6, ET0.wdc, coords)$corr, modified.ttest(ET0.cmip6, ET0.wdc, coords)$p.value)))

df.combined <- rbind(
  cbind(Metric = c("r","p-value","r2","Dutilleul","p-value"), df.era5.wdc),
  cbind(Metric = c("r","p-value","r2","Dutilleul","p-value"), df.era5.cmip6),
  cbind(Metric = c("r","p-value","r2","Dutilleul","p-value"), df.wdc.cmip6)
)

knitr::kable(df.combined, 
             digits = 2, caption = "Pearson correlation coefficients and r2 between ERA5, Worldclim and CMIP6") %>%
  kableExtra::kable_styling(bootstrap_options = "bordered") %>%
  kableExtra::group_rows(group_label = "ERA5 vs WDC", start_row = 1, end_row = 5) %>%
  kableExtra::group_rows("ERA5 vs CMIP6", 6,10) %>%
  kableExtra::group_rows("WDC vs CMIP6", 11,15)

Pearson correlation coefficients and r2 between ERA5, Worldclim and CMIP6

Metric	precipitation	temperature	wind	Rn	ea	ET0
ERA5 vs WDC
r	0.88	1.00	0.77	0.89	1.00	0.97
p-value	0.00	0.00	0.00	0.00	0.00	0.00
r2	0.77	1.00	0.59	0.79	0.99	0.94
Dutilleul	0.88	1.00	0.77	0.89	1.00	0.97
p-value	0.00	0.00	0.00	0.00	0.00	0.00
ERA5 vs CMIP6
r	0.89	1.00	0.88	0.99	0.99	0.98
p-value	0.00	0.00	0.00	0.00	0.00	0.00
r2	0.79	1.00	0.77	0.97	0.98	0.96
Dutilleul	0.89	1.00	0.88	0.99	0.99	0.98
p-value	0.00	0.00	0.00	0.00	0.00	0.00
WDC vs CMIP6
r	0.89	1.00	0.78	0.91	0.99	0.96
p-value	0.00	0.00	0.00	0.00	0.00	0.00
r2	0.79	0.99	0.61	0.82	0.97	0.91
Dutilleul	0.89	1.00	0.78	0.91	0.99	0.96
p-value	0.00	0.00	0.00	0.00	0.00	0.00

5.2.2 ERA5 compared to WDC

Rn is computed differently for ERA5 and WDC, hence the discrepancy for low values.

ggarrange(

ggplot(calib.noant,aes(x = tas.wdc))+geom_point(aes( y = tas.era5), alpha =0.1, shape = 16, size = 0.5, col = colorvec[1])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, ERA5")+
  lims(x = c(-30, 31), y = c(-30,31))+
   annotate(geom = "text", x =  min(c(-30, 31)) + dist(c(-30, 31))*0.15, y =  max(c(-30, 31)) - dist(c(-30, 31))*0.15, 
            label = paste("r = ", round(cor(calib.noant$tas.wdc, calib.noant$tas.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5), alpha =0.1, shape = 16, size = 0.5, col = colorvec[1])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  lims(x = c(-1, 6000), y = c(-1, 6000))+
   annotate(geom = "text", x =  min(c(-1, 6000)) + dist(c(-1, 6000))*0.15, y =  max(c(-1, 6000)) - dist(c(-1, 6000))*0.15, 
            label = paste("r = ", round(cor(calib.noant$pr.wdc, calib.noant$pr.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5), alpha =0.1, shape = 16, size = 0.5, col = colorvec[1])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5")+
  lims(x = c(0,10), y = c(0,10))+
   annotate(geom = "text", x =  min(c(0,10)) + dist(c(0,10))*0.15, y =  max(c(0,10)) - dist(c(0,10))*0.15,
            label = paste("r = ", round(cor(calib.noant$sfcWind.wdc, calib.noant$sfcWind.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5), alpha = 0.1, shape = 16, size = 0.5, col = colorvec[1])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   lims(x = c(-8,20), y = c(-8,20))+
   annotate(geom = "text", x =  min(c(-8,20)) + dist(c(-8,20))*0.15, y =  max(c(-8,20)) - dist(c(-8,20))*0.15,
            label = paste("r = ", round(cor(calib.noant$Rn.wdc, calib.noant$Rn.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[1])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC", y = "ea, kPa, ERA5", col = "Dataset")+
   lims(x = c(0,3), y = c(0,3))+
   annotate(geom = "text",  x =  min(c(0,3)) + dist(c(0,3))*0.15, y =  max(c(0,3)) - dist(c(0,3))*0.15,
            label = paste("r = ", round(cor(calib.noant$ea.wdc, calib.noant$ea.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[1])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
lims(x = c(-261,4200), y = c(-261,4200))+
   annotate(geom = "text",  x =  min(c(-261,4200)) + dist(c(-261,4200))*0.15, y =  max(c(-261,4200)) - dist(c(-261,4200))*0.15,
            label = paste("r = ", round(cor(calib.noant$ET0.wdc, calib.noant$ET0.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

ggplot(calib.noant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none")


ggarrange(

ggplot(calib.noant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

5.2.3 CMIP6 and WDC

ggarrange(

ggplot(calib.noant,aes(x = tas.wdc))+geom_point(aes( y = tas.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[4])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, CMIP6")+
  lims(x = c(-30, 31), y = c(-30,31))+
   annotate(geom = "text", x =  min(c(-30, 31)) + dist(c(-30, 31))*0.15, y =  max(c(-30, 31)) - dist(c(-30, 31))*0.15, 
            label = paste("r = ", round(cor(calib.noant$tas.wdc, calib.noant$tas.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = pr.wdc))+geom_point(aes( y = pr.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[4])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, CMIP6")+
  scale_color_manual(values = c(colorvec[4], colorvec[4]))+
  lims(x = c(-1, 6000), y = c(-1, 6000))+
   annotate(geom = "text", x =  min(c(-1, 6000)) + dist(c(-1, 6000))*0.15, y =  max(c(-1, 6000)) - dist(c(-1, 6000))*0.15, 
            label = paste("r = ", round(cor(calib.noant$pr.wdc, calib.noant$pr.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[4])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, CMIP6")+
  lims(x = c(0,10), y = c(0,10))+
   annotate(geom = "text", x =  min(c(0,10)) + dist(c(0,10))*0.15, y =  max(c(0,10)) - dist(c(0,10))*0.15,
            label = paste("r = ", round(cor(calib.noant$sfcWind.wdc, calib.noant$sfcWind.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.cmip6), alpha = 0.1, shape = 16, size = 0.5, col = colorvec[4])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, CMIP6", col = "Dataset")+
   lims(x = c(-8,20), y = c(-8,20))+
   annotate(geom = "text", x =  min(c(-8,20)) + dist(c(-8,20))*0.15, y =  max(c(-8,20)) - dist(c(-8,20))*0.15,
            label = paste("r = ", round(cor(calib.noant$Rn.wdc, calib.noant$Rn.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ea.wdc))+geom_point(aes(y = ea.cmip6, col = "cmip6"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[4])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, CMIP6", col = "Dataset")+
   lims(x = c(0,3), y = c(0,3))+
   annotate(geom = "text",  x =  min(c(0,3)) + dist(c(0,3))*0.15, y =  max(c(0,3)) - dist(c(0,3))*0.15,
            label = paste("r = ", round(cor(calib.noant$ea.wdc, calib.noant$ea.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.cmip6, col = "cmip6"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[4])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, CMIP6", col = "Dataset")+
lims(x = c(-261,4200), y = c(-261,4200))+
   annotate(geom = "text",  x =  min(c(-261,4200)) + dist(c(-261,4200))*0.15, y =  max(c(-261,4200)) - dist(c(-261,4200))*0.15,
            label = paste("r = ", round(cor(calib.noant$ET0.wdc, calib.noant$ET0.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

5.2.4 CMIP6 and ERA5

ggarrange(
  
ggplot(calib.noant,aes(x = tas.era5))+geom_point(aes( y = tas.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[8])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, ERA5", y = "Temperature in °C, CMIP6")+
  lims(x = c(-30, 31), y = c(-30,31))+
   annotate(geom = "text", x =  min(c(-30, 31)) + dist(c(-30, 31))*0.15, y =  max(c(-30, 31)) - dist(c(-30, 31))*0.15, 
            label = paste("r = ", round(cor(calib.noant$tas.era5, calib.noant$tas.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = pr.era5))+geom_point(aes( y = pr.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[8])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, ERA5", y = "Annual precipitation, mm, CMIP6")+
  scale_color_manual(values = c(colorvec[4], colorvec[4]))+
  lims(x = c(-1, 6000), y = c(-1, 6000))+
   annotate(geom = "text", x =  min(c(-1, 6000)) + dist(c(-1, 6000))*0.15, y =  max(c(-1, 6000)) - dist(c(-1, 6000))*0.15, 
            label = paste("r = ", round(cor(calib.noant$pr.era5, calib.noant$pr.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = sfcWind.era5))+geom_point(aes( y = sfcWind.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[8])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s, CMIP6")+
  lims(x = c(0,10), y = c(0,10))+
   annotate(geom = "text", x =  min(c(0,10)) + dist(c(0,10))*0.15, y =  max(c(0,10)) - dist(c(0,10))*0.15,
            label = paste("r = ", round(cor(calib.noant$sfcWind.era5, calib.noant$sfcWind.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = Rn.era5))+geom_point(aes( y = Rn.cmip6), alpha = 0.1, shape = 16, size = 0.5, col = colorvec[8])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, ERA5, °C", y = "Rn - G, MJ/m2/day, CMIP6", col = "Dataset")+
   lims(x = c(-8,20), y = c(-8,20))+
   annotate(geom = "text", x =  min(c(-8,20)) + dist(c(-8,20))*0.15, y =  max(c(-8,20)) - dist(c(-8,20))*0.15,
            label = paste("r = ", round(cor(calib.noant$Rn.era5, calib.noant$Rn.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ea.era5))+geom_point(aes(y = ea.cmip6, col = "cmip6"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[8])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, ERA5, °C", y = "ea, kPa, CMIP6", col = "Dataset")+
   lims(x = c(0,3), y = c(0,3))+
   annotate(geom = "text",  x =  min(c(0,3)) + dist(c(0,3))*0.15, y =  max(c(0,3)) - dist(c(0,3))*0.15,
            label = paste("r = ", round(cor(calib.noant$ea.era5, calib.noant$ea.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ET0.era5))+geom_point(aes(y = ET0.cmip6, col = "cmip6"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[8])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, ERA5", y = "ET0, mm/day, CMIP6", col = "Dataset")+
lims(x = c(-261,4200), y = c(-261,4200))+
   annotate(geom = "text",  x =  min(c(-261,4200)) + dist(c(-261,4200))*0.15, y =  max(c(-261,4200)) - dist(c(-261,4200))*0.15,
            label = paste("r = ", round(cor(calib.noant$ET0.era5, calib.noant$ET0.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2

)

5.3 Scatterplots only Antarctica

calib.ant <- calib.wide %>% filter(lat < -60)

5.3.1 ERA5 and WDC

Rn is computed differently for ERA5 and WDC, hence the discrepancy for low values.

ggarrange(

ggplot(calib.ant,aes(x = tas.wdc))+geom_point(aes( y = tas.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

ggplot(calib.ant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none")


ggarrange(

ggplot(calib.ant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

5.3.2 CMIP6 and WDC

ggarrange(

ggplot(calib.ant,aes(x = tas.wdc))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = pr.wdc))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = sfcWind.wdc))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = Rn.wdc))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ea.wdc))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ET0.wdc))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2
)

5.3.3 CMIP6 and ERA5

ggarrange(
  
ggplot(calib.ant,aes(x = tas.era5))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, ERA5", y = "Temperature in °C, CMIP6", col = "Dataset")+
  scale_color_manual(values = colorvec[8])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = pr.era5))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, ERA5", y = "Annual precipitation, mm, CMIP6", col = "Dataset")+
  scale_color_manual(values = c(colorvec[8]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = sfcWind.era5))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s, CMIP6", col = "Dataset")+
  scale_color_manual(values = c(colorvec[8]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = Rn.era5))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
   stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, ERA5, °C", y = "Rn - G, W/MJ/m2/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ea.era5))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, ERA5", y = "ea, kPa, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ET0.era5))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, ERA5", y = "ET0, mm/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2)

5.4 Aridity category difference

5.4.1 ERA5 and WDC

calib.wide$same.cat.wdc <- with(calib.wide, ifelse(cat.AI.era5 == cat.AI.wdc, "y","n"))
calib.wide$change.cat.wdc <- with(calib.wide, paste(cat.AI.era5, cat.AI.wdc, sep = " to "))
table(calib.wide$change.cat.wdc)
## 
##                 Arid to Arid                Arid to Humid 
##                         1359                            1 
##            Arid to Hyperarid            Arid to Semi-arid 
##                           18                          155 
##                 Cold to Arid                 Cold to Cold 
##                           19                        10101 
##         Cold to Dry subhumid                Cold to Humid 
##                            9                          290 
##            Cold to Semi-arid         Dry subhumid to Arid 
##                           20                           23 
##         Dry subhumid to Cold Dry subhumid to Dry subhumid 
##                            5                          205 
##        Dry subhumid to Humid    Dry subhumid to Semi-arid 
##                          209                          271 
##                Humid to Arid                Humid to Cold 
##                           30                          520 
##        Humid to Dry subhumid               Humid to Humid 
##                          347                         5436 
##           Humid to Semi-arid            Hyperarid to Arid 
##                          210                          230 
##            Hyperarid to Cold       Hyperarid to Hyperarid 
##                            2                          726 
##                   NA to Cold            Semi-arid to Arid 
##                          570                          181 
##            Semi-arid to Cold    Semi-arid to Dry subhumid 
##                            3                          173 
##           Semi-arid to Humid       Semi-arid to Semi-arid 
##                           71                         1077
calib.wide$dryer.wdc <- with(calib.wide, ifelse(change.cat.wdc == "Arid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Arid to Semi-arid", "dryer",
                                ifelse(change.cat.wdc == " Arid to Dry subhumid", "dryer", 
                                ifelse(change.cat.wdc == "Arid to Hyperarid", "wetter",
                                ifelse(change.cat.wdc == "Arid to Cold", "dryer",       
                                ifelse(change.cat.wdc == "Cold to Humid", "dryer",
                                ifelse(change.cat.wdc == "Cold to Hyperarid", "wetter",
                                ifelse(change.cat.wdc == "Dry subhumid to Cold", "dryer",
                                ifelse(change.cat.wdc == "Dry subhumid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Dry subhumid to Arid", "wetter",
                                ifelse(change.cat.wdc == "Dry subhumid to Semi-arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to Semi-arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to Cold","dryer",
                                ifelse(change.cat.wdc == "Humid to Dry subhumid","wetter",
                                ifelse(change.cat.wdc == "Semi-arid to Dry subhumid", "dryer",
                                ifelse(change.cat.wdc == "Semi-arid to Arid", "wetter",
                                ifelse(change.cat.wdc == "Semi-arid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Semi-arid to Cold", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Arid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Semi-arid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Dry subhumid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Humid", "dryer",
                                NA))))))))))))))))))))))))

ggarrange(
ggplot(filter(calib.wide, same.cat.wdc == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.wdc))+
  scale_fill_manual(values = c( "#ef8a62", "#67a9cf"), na.translate = F)+
  borders()+
  labs(fill = "ERA5 compared to WDC", title = "Change of category")+
  ylim(-55,90)+
  theme_void(base_size = 15)+
  theme(legend.position = "bottom"),

ggplot() + geom_raster(data = calib.wide, aes(x=lon, y = lat,  fill = AI.era5-AI.wdc))+
  binned_scale(aesthetics = "fill", breaks = c(-1, -0.1, -0.01, 0, 0.01, 0.1, 1), 
               palette = function(x) rev(c("#08519c", "#2166ac", "#67a9cf", "aliceblue","bisque", "#ef8a62", "#b2182b", "#67001f")),
               guide = guide_legend(label.theme = element_text(angle = 0))
               )+
  labs(title = "Changes in AI", fill = "ERA5 compared to WDC")+
  theme_void(base_size = 15)+ylim(-55,90)+
  theme(legend.position = "bottom"),
ncol = 2)

ggplot() + 
  geom_raster(data = subset(calib.wide, pr.era5 - pr.wdc < 277 & pr.era5 - pr.wdc > - 98), aes(x=lon, y = lat,  fill = pr.era5-pr.wdc))+
  geom_raster(data = subset(calib.wide, pr.era5 - pr.wdc > 277), aes(x=lon, y = lat), fill = "#44FF44")+
   geom_raster(data = subset(calib.wide, pr.era5 - pr.wdc < - 98), aes(x=lon, y = lat), fill = "#FF4444")+
  labs(title = "Most extreme differences in precipitation", subtitle = "In red, the 10% most extreme gricells in which Worldclim is wetter than ERA5; \nin green, the 10% most extreme gridcells in which ERA5 is wetter than Worldclim", fill = "precipitation in ERA5 - precipitations in WDC,\nmm/year")+
  theme_void(base_size = 15)+
  theme(legend.position = "bottom")

5.4.2 CMIP6 and WDC

calib.wide$same.cat.wdc <- with(calib.wide, ifelse(cat.AI.cmip6 == cat.AI.wdc, "y","n"))
calib.wide$change.cat.wdc <- with(calib.wide, paste(cat.AI.cmip6, cat.AI.wdc, sep = " to "))
table(calib.wide$change.cat.wdc)
## 
##                 Arid to Arid         Arid to Dry subhumid 
##                         1099                            9 
##            Arid to Hyperarid            Arid to Semi-arid 
##                           83                          150 
##                 Cold to Arid                 Cold to Cold 
##                           25                        10491 
##         Cold to Dry subhumid                Cold to Humid 
##                            5                          588 
##            Cold to Semi-arid         Dry subhumid to Arid 
##                           23                           36 
## Dry subhumid to Dry subhumid        Dry subhumid to Humid 
##                          147                          196 
##    Dry subhumid to Hyperarid    Dry subhumid to Semi-arid 
##                            1                          350 
##                Humid to Arid                Humid to Cold 
##                           77                          138 
##        Humid to Dry subhumid               Humid to Humid 
##                          427                         5112 
##           Humid to Hyperarid           Humid to Semi-arid 
##                            1                          435 
##            Hyperarid to Arid            Hyperarid to Cold 
##                           85                            2 
##       Hyperarid to Hyperarid                   NA to Cold 
##                          652                          570 
##            Semi-arid to Arid    Semi-arid to Dry subhumid 
##                          520                          146 
##           Semi-arid to Humid       Semi-arid to Hyperarid 
##                          111                            7 
##       Semi-arid to Semi-arid 
##                          775
calib.wide$dryer.wdc <- with(calib.wide, ifelse(change.cat.wdc == "Arid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Arid to Semi-arid", "dryer",
                                ifelse(change.cat.wdc == " Arid to Dry subhumid", "dryer", 
                                ifelse(change.cat.wdc == "Arid to Hyperarid", "wetter",
                                ifelse(change.cat.wdc == "Arid to Cold", "dryer",       
                                ifelse(change.cat.wdc == "Cold to Humid", "dryer",
                                ifelse(change.cat.wdc == "Cold to Hyperarid", "wetter",
                                ifelse(change.cat.wdc == "Dry subhumid to Cold", "dryer",
                                ifelse(change.cat.wdc == "Dry subhumid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Dry subhumid to Arid", "wetter",
                                ifelse(change.cat.wdc == "Dry subhumid to Semi-arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to Semi-arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to Cold","dryer",
                                ifelse(change.cat.wdc == "Humid to Dry subhumid","wetter",
                                ifelse(change.cat.wdc == "Semi-arid to Dry subhumid", "dryer",
                                ifelse(change.cat.wdc == "Semi-arid to Arid", "wetter",
                                ifelse(change.cat.wdc == "Semi-arid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Semi-arid to Cold", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Arid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Semi-arid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Dry subhumid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Humid", "dryer",
                                NA))))))))))))))))))))))))

ggarrange(
ggplot(filter(calib.wide, same.cat.wdc == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.wdc))+
  scale_fill_manual(values = c( "#ef8a62", "#67a9cf"), na.translate = F)+
  borders()+
  labs(fill = "CMIP6 compared to WDC", "Change of category")+
  ylim(-55,90)+
  theme_void(base_size = 15)+
  theme(legend.position = "bottom"),

ggplot() + geom_raster(data = calib.wide, aes(x=lon, y = lat,  fill = AI.cmip6-AI.wdc))+
   borders()+
  binned_scale(aesthetics = "fill", breaks = c(-1, -0.1, -0.01, 0, 0.01, 0.1, 1), palette = function(x) rev(c("#08519c", "#2166ac", "#67a9cf", "aliceblue","bisque", "#ef8a62", "#b2182b", "#67001f")), 
                                                                                                              guide = guide_legend(label.theme = element_text(angle = 0)))+
  labs(title = "Change in AI", fill = "CMIP6 compared to WDC")+
  theme_void(base_size = 15)+ylim(-55,90)+
  theme(legend.position = "bottom"),

ncol = 2)

5.5 CMIP6 and ERA5

calib.wide$same.cat.era5 <- with(calib.wide, ifelse(cat.AI.cmip6 == cat.AI.era5, "y","n"))
calib.wide$change.cat.era5 <- with(calib.wide, paste(cat.AI.cmip6, cat.AI.era5, sep = " to "))
table(calib.wide$change.cat.era5)
## 
##                 Arid to Arid         Arid to Dry subhumid 
##                          989                            6 
##                Arid to Humid            Arid to Hyperarid 
##                            2                          230 
##            Arid to Semi-arid                 Cold to Cold 
##                          114                        10377 
##         Cold to Dry subhumid                Cold to Humid 
##                           10                          745 
##         Dry subhumid to Arid Dry subhumid to Dry subhumid 
##                           14                          169 
##        Dry subhumid to Humid    Dry subhumid to Semi-arid 
##                          219                          328 
##                Humid to Arid                Humid to Cold 
##                            5                           62 
##        Humid to Dry subhumid               Humid to Humid 
##                          382                         5456 
##           Humid to Semi-arid            Hyperarid to Arid 
##                          285                           14 
##       Hyperarid to Hyperarid                     NA to NA 
##                          725                          570 
##            Semi-arid to Arid    Semi-arid to Dry subhumid 
##                          511                          146 
##           Semi-arid to Humid       Semi-arid to Hyperarid 
##                          121                            3 
##       Semi-arid to Semi-arid 
##                          778
calib.wide$dryer.era5 <- with(calib.wide, ifelse(change.cat.era5 == "Arid to Humid", "dryer",
                                ifelse(change.cat.era5 == "Arid to Semi-arid", "dryer",
                                ifelse(change.cat.era5 == " Arid to Dry subhumid", "dryer", 
                                ifelse(change.cat.era5 == "Arid to Hyperarid", "wetter",
                                ifelse(change.cat.era5 == "Cold to Humid", "dryer",
                                ifelse(change.cat.era5 == "Dry subhumid to Cold", "dryer",
                                ifelse(change.cat.era5 == "Dry subhumid to Humid", "dryer",
                                ifelse(change.cat.era5 == "Dry subhumid to Arid", "wetter",
                                ifelse(change.cat.era5 == "Dry subhumid to Semi-arid", "wetter",
                                ifelse(change.cat.era5 == "Humid to arid", "wetter",
                                ifelse(change.cat.era5 == "Humid to Semi-arid", "wetter",
                                ifelse(change.cat.era5 == "Humid to Cold","dryer",
                                ifelse(change.cat.era5 == "Humid to Dry subhumid","wetter",
                                ifelse(change.cat.era5 == "Semi-arid to Dry subhumid", "dryer",
                                ifelse(change.cat.era5 == "Semi-arid to Arid", "wetter",
                                ifelse(change.cat.era5 == "Semi-arid to Humid", "dryer",
                                ifelse(change.cat.era5 == "Hyperarid to Arid", "dryer",
                                ifelse(change.cat.era5 == "Hyperarid to Semi-arid", "dryer",
                                NA)))))))))))))))))))

ggarrange(
ggplot(filter(calib.wide, same.cat.era5 == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.era5))+
  scale_fill_manual(values = c( "#ef8a62", "#67a9cf"), na.translate = F)+
  borders()+
  labs(fill = "CMIP6 compared to ERA5", "Change of aridity category")+
  ylim(-55,90)+
  theme_void()+
  theme(legend.position = "bottom"),


ggplot() + geom_raster(data = calib.wide, aes(x=lon, y = lat,  fill = AI.cmip6-AI.era5))+
   borders()+
  binned_scale(aesthetics = "fill", breaks =c(-1, -0.1, -0.01, 0, 0.01, 0.1, 1), palette = function(x) rev(c("#08519c", "#2166ac", "#67a9cf", "aliceblue","bisque", "#ef8a62", "#b2182b", "#67001f")))+
  labs(title = "Change of AI", fill = "CMIP6 compared to ERA5")+
  theme_void()+ylim(-55,90)+
  theme(legend.position = "bottom"),
ncol = 2)

6 NDVI

# Download Gimms NDVI
library(gimms)
library(gganimate)

land_mask <- raster("Masks/land_sea_mask_1degree.nc4")
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land

#ndvi.ref <- downloadGimms(x= 1981,y= 2000,dsn = "NDVI")

list.ndvi <- list.files(path = "NDVI", pattern = "ndvi.*\\.nc4$", full.names = T)[1:2]

ndvi.rast <- rasterizeGimms(list.ndvi, keep = 0) %>% mean(na.rm = T) %>% projectRaster(land_mask) # rasterize Gimms does the quality control and removes NAs
ndvi.df <- ndvi.rast %>% as.data.frame(xy = T) %>% setNames(c("lon","lat","ndvi")) %>% merge(land_mask.df, by = c("lon","lat")) %>% filter(lm != 0)

write.table(ndvi.df, "NDVI/ndvi.df.txt")

ndvi.df <- read.table("NDVI/ndvi.df.txt") 

ggplot(subset(ndvi.df, lm == 1))+geom_raster(aes(x=lon, y = lat, fill = ndvi))+
  scale_fill_continuous_sequential(palette = "Greens 3", rev = T)+
  ylim(-55,90)+labs(fill = "NDVI")+
  theme_void()+theme(legend.position = "bottom")

calibn <- merge(calib, ndvi.df, by = c("lon","lat"))

ggplot(subset(calibn, !is.na(cat.AI)))+geom_boxplot(aes(x=cat.AI, y = ndvi, col = cat.AI))+facet_grid(rows = vars(source))+
  scale_color_manual(values = col.cat)+
  labs(x = "", y = "NDVI")+
  theme_minimal()+theme(legend.position = "none")

aov(ndvi~cat.AI, data = group_by(calibn, source)) %>% summary()
##                Df Sum Sq Mean Sq F value Pr(>F)    
## cat.AI          5   1026  205.17    8708 <2e-16 ***
## Residuals   43200   1018    0.02                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 23577 observations effacées parce que manquantes

range.ndvi <- calibn %>% filter(!is.na(cat.AI)) %>% group_by(source, cat.AI) %>% 
  summarise(range.ndvi = paste(round(quantile(ndvi, na.rm = T)[c(2,4)],2), collapse = " to ")) %>%
  cast(source~cat.AI)

cor.ndvi <- calibn %>% filter(!is.na(cat.AI)) %>% group_by(source, cat.AI) %>% 
  summarise(pearson.r = round(cor(ndvi, AI, use = "pairwise.complete"),2)) %>%
  cast(source~cat.AI)

knitr::kable(range.ndvi, caption = "NDVI range by source and aridity category")

NDVI range by source and aridity category

source	Cold	Hyperarid	Arid	Semi-arid	Dry subhumid	Humid
CMIP6	0.3 to 0.47	0.09 to 0.11	0.1 to 0.19	0.17 to 0.36	0.24 to 0.49	0.4 to 0.69
ERA5	0.27 to 0.47	0.09 to 0.11	0.11 to 0.21	0.19 to 0.37	0.24 to 0.48	0.42 to 0.68
Worldclim	0.29 to 0.47	0.09 to 0.11	0.1 to 0.19	0.18 to 0.34	0.28 to 0.44	0.45 to 0.7

knitr::kable(cor.ndvi, caption = "Correlation between NDVI and AI by source and aridity category")

Correlation between NDVI and AI by source and aridity category

source	Cold	Hyperarid	Arid	Semi-arid	Dry subhumid	Humid
CMIP6	-0.08	0.15	0.36	0.38	-0.01	0.31
ERA5	-0.16	0.08	0.47	0.30	0.08	0.38
Worldclim	0.05	0.05	0.49	0.39	0.14	0.52

knitr::knit_exit()

6.0.1 Maps CMIP6 ERA5 // WDC

In red: when CMIP6 or ERA5 are warmer than Wordlclim

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.era5-tas.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: ERA5 - WDC")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.cmip6-tas.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: CMIP6 - WDC")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when CMIP6 or ERA5 are wetter than Worldclim

cowplot::plot_grid(
  
  ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.era5-pr.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: ERA5 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.cmip6-pr.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: CMIP6 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when wind speed is higher for CMIP6 or ERA5 compared to Worldclim

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.era5-sfcWind.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: ERA5 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.cmip6-sfcWind.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: CMIP6 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)

6.0.2 Maps CMIP6 WDC // ERA5

In red: when CMIP6 or WDC are warmer than ERA5

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.wdc-tas.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: WDC - ERA5")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.cmip6-tas.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: CMIP6 - ERA5")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when CMIP6 or WDC are wetter than ERA5

cowplot::plot_grid(
  
  ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.wdc-pr.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: WDC - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.cmip6-pr.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: CMIP6 - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when wind speed is higher for CMIP6 or WDC compared to ERA5

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.wdc-sfcWind.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: WDC - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.cmip6-sfcWind.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: CMIP6 - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)